#include <assert.h>
#include <petscksp.h>
#include <petscmat.h>

int main(int argc,
	 char** argv)
{
    PetscInitialize(&argc, &argv, (char *)0, (char *)0);
    PetscPrintf(PETSC_COMM_WORLD, "%s", "... Start\n");
    
    Vec x; /* approx solution, RHS, exact solution */
    Vec b;
    Vec u;
    Mat matA; /* linear system matrix */
    KSP ksp; /* linear solver context */
    PC pc;  /* preconditioner context */
    PetscReal      norm;         /* norm of solution error */
    PetscInt       its;        //n = 10, col[3]
    PetscErrorCode ierr;

    /**
     * Reads CRS format into PETSC SEQAIJ format
     */
    PetscInt nr;
    PetscInt nv;

    /* read case0.dat to matrix dimensions */
    FILE* finCaseDat0;
    finCaseDat0 = fopen("case.dat0", "r");
    fscanf(finCaseDat0, "%d", &nr);
    fscanf(finCaseDat0, "%d", &nv);
    assert(nr == 9);
    assert(nv == 5);
    close(finCaseDat0);
    PetscInt n = nv*nr;
    
    CHKERRQ(ierr = MatCreate(PETSC_COMM_WORLD, &matA));
    
    CHKERRQ(ierr = MatSetSizes(matA,
			       PETSC_DECIDE,
			       PETSC_DECIDE,
			       n,
			       n));
    
    CHKERRQ(ierr = MatSetFromOptions(matA));

    /* Create right-hand-side vector b */
    CHKERRQ(ierr = VecCreate(PETSC_COMM_WORLD, &b));

    CHKERRQ(ierr = VecSetSizes(b, PETSC_DECIDE, n));

    CHKERRQ(ierr = VecSetFromOptions(b));

    CHKERRQ(ierr = VecDuplicate(b, &x));

    FILE* finRHS = fopen("rhs_value.dat", "r");
    PetscInt rhsi;
    for (rhsi = 0; rhsi < n; ++rhsi)
    {
	assert(!feof(finRHS));
	PetscScalar v;	
	fscanf(finRHS, "%le", &v);
	printf("%le", v);
	CHKERRQ(ierr = VecSetValue(b, rhsi, v, INSERT_VALUES));
    }
    close(finRHS);

    CHKERRQ(ierr = VecDuplicate(b,&x));
    CHKERRQ(ierr = VecDuplicate(b,&u));

    CHKERRQ(ierr = VecAssemblyBegin(b));
    CHKERRQ(ierr = VecAssemblyEnd(b));

    CHKERRQ(ierr = VecAssemblyBegin(u));
    CHKERRQ(ierr = VecAssemblyEnd(u));
    
    CHKERRQ(ierr = VecAssemblyBegin(x));
    CHKERRQ(ierr = VecAssemblyEnd(x));

    CHKERRQ(ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF));
    
    
    FILE* finV = fopen("matrix_coef.dat", "r");
    FILE* finR = fopen("row_ptr.dat", "r"); 
    FILE* finC = fopen("col_ind.dat", "r");

    PetscInt row = 0;   /* row counter */
    PetscInt rs;  /* start index in row_ptr.dat */
    PetscInt re;  /* end index in row_ptr.dat */
    fscanf(finR, "%d", &rs);
    fscanf(finR, "%d", &re);
    while(!feof(finR))
    {
	PetscInt range = re - rs;
	//printf("%d", range);
	PetscInt i;
	for (i = 0; i < range; ++i)
	{
	    assert(!feof(finC));
	    PetscInt col;
	    fscanf(finC, "%d", &col);
	    
	    PetscInt j;
	    for (j = 0; j < nv; ++j)
	    {
		PetscInt cols[nv];
		PetscScalar vals[nv];
		PetscInt k;
		for (k = 0; k < nv; k++)
		{
		    assert(!feof(finV));
		    PetscReal v;
		    fscanf(finV, "%le", &v);
		    //printf("%le", &v);
		    cols[k] = col*nv+k;
		    vals[k] = v;
		}
		PetscInt itmp = row*nv+j;
		ierr = MatSetValues(matA, 1, &itmp, 5, cols, vals, INSERT_VALUES); 
	    }
	}
	
	/* update to another row */
	rs = re;
	assert(!feof(finR));
	fscanf(finR, "%d", &re);
	++row;
    }

    close(finV);
    close(finR);
    close(finC);

    CHKERRQ(ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY));

    CHKERRQ(ierr = MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY));
  
    CHKERRQ(ierr = MatView(matA, PETSC_VIEWER_STDOUT_SELF));    
   
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "...done\n"));

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
       Create the linear solver and set various options
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    /* 
       Create linear solver context
    */
    CHKERRQ(ierr = KSPCreate(PETSC_COMM_WORLD,
			 &ksp));

    /* 
       Set operators. Here the matrix that defines the linear system
       also serves as the preconditioning matrix.
    */
    CHKERRQ(ierr = KSPSetOperators(ksp,
				   matA,
				   matA,
				   DIFFERENT_NONZERO_PATTERN));
    
    /* 
       Set linear solver defaults for this problem (optional).
       - By extracting the KSP and PC contexts from the KSP context,
       we can then directly call any KSP and PC routines to set
       various options.
       - The following four statements are optional; all of these
       parameters could alternatively be specified at runtime via
       KSPSetFromOptions();
    */
    CHKERRQ(ierr = KSPGetPC(ksp, &pc));
   
    CHKERRQ(ierr = PCSetType(pc, PCLU));
    
    //    ierr = PCFactorSetMatSolverPackage(pc,
    //				       "mumps");
    //CHKERRQ(ierr);

    CHKERRQ(ierr = KSPSetTolerances(ksp,
				    1.e-7,
				    PETSC_DEFAULT,
				    PETSC_DEFAULT,
				    PETSC_DEFAULT));
    
    
    /* 
       Set runtime options, e.g.,
       -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
       These options will override those specified above as long as
       KSPSetFromOptions() is called _after_ any other customization
       routines.
    */
    CHKERRQ(ierr = KSPSetFromOptions(ksp));
    
    /* KSP Set Type to KSPGMRES */
    CHKERRQ(ierr = KSPSetType(ksp, KSPGMRES));
    
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
       Solve the linear system
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    /* 
       Solve linear system
    */
    CHKERRQ(ierr = KSPSolve(ksp,b,x));
     

/*
    VecAssemblyBegin(x);
    CHKERRQ(ierr);
    VecAssemblyEnd(x);
    CHKERRQ(ierr);
*/

    PetscInt one = 1;
    CHKERRQ(ierr = VecSet(u,one));
    

/*    
    VecAssemblyBegin(u);
    CHKERRQ(ierr);
    VecAssemblyEnd(u);
    CHKERRQ(ierr);
*/    
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));
    CHKERRQ(ierr = VecView(x, PETSC_VIEWER_STDOUT_SELF));
    
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));
    CHKERRQ(ierr = VecView(u, PETSC_VIEWER_STDOUT_SELF));

    CHKERRQ(ierr = PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));

    /* 
       View solver info; we could instead use the option -ksp_view to
       print this info to the screen at the conclusion of KSPSolve().
    */
    CHKERRQ(ierr = KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));

    /* MatDestroy(matA) */
    /* VecDestroy(b) */
    /* VecDestroy(x) */

    CHKERRQ(VecNorm(u,NORM_2,&norm));
    CHKERRQ(KSPGetIterationNumber(ksp,&its));
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,
			"Norm of error %A, Iterations %D\n",
			norm,its));



    PetscFinalize();

    return 0;
}
